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Abstract. We develop a perturbative approach to redshift space distortions (RSD) using the phase 
space distribution function approach and apply it to the dark matter redshift space power spectrum 
and its moments. RSD can be written as a sum over density weighted velocity moments correlators, 
with the lowest order being density, momentum density and stress energy density. We use standard 
and extended perturbation theory (PT) to determine their auto and cross correlators, comparing 
them to N-body simulations. We show which of the terms can be modeled well with the standard 
PT and which need additional terms that include higher order corrections which cannot be modeled 
in PT. Most of these additional terms are related to the small scale velocity dispersion effects, the 
so called finger of god (FoG) effects, which affect some, but not all, of the terms in this expansion, 
and which can be approximately modeled using a simple physically motivated ansatz such as the halo 
model. We point out that there are several velocity dispersions that enter into the detailed RSD 
analysis with very different amplitudes, which can be approximately predicted by the halo model. In 
contrast to previous models our approach systematically includes all of the terms at a given order in 
PT and provides a physical interpretation for the small scale dispersion values. We investigate RSD 
power spectrum as a function of fi, the cosine of the angle between the Fourier mode and line of sight, 
focusing on the lowest order powers of fi and multipole moments which dominate the observable RSD 
power spectrum. Overall we find considerable success in modeling many, but not all, of the terms 
in this expansion. This is similar to the situation in real space, but predicting power spectrum in 
redshift space is more difficult because of the explicit influence of small scale dispersion type effects 
in RSD, which extend to very large scales. 
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1 Introduction 

Galaxy clustering surveys are one of the most important venues to extract cosmological information 
today. This is because by measuring the 3 dimensional distribution of galaxies we can in principle 
relate it to the 3 dimensional distribution of the underlying dark matter. The dark matter distribution 
is sensitive to many of the cosmological parameters. The growth of dark matter structures in time 
also provides important constraints on the models, such as the nature and amount of dark energy. 

Since galaxies are not perfect tracers of dark matter, their clustering is biased relative to the dark 
matter. This means that galaxy surveys cannot determine the rate of growth of structure unless this 
biasing is determined. Fortunately, galaxy redshift surveys provide additional information, because 
the observed redshift is a sum of the radial distance to the galaxy and its peculiar velocity (Dopplcr 
shift). Galaxies are expected to follow the same gravitational potential as the dark matter and thus 
they are expected to have the same velocity (in a large-scale average at least). This leads to a 
clustering strength that depends on the angle between the galaxy pairs and the line of sight, which 
is referred to as redshift space distortions (RSD). In linear theory it can be easily related to the dark 
matter clustering [1, 2]. These distortions thus make the galaxy clustering in redshift space more 
complex, but at the same time provide an opportunity to extract important information on the dark 
matter clustering directly from the redshift surveys. To what extent this is possible is a matter of 
considerable debate: there are significant nonlinear effects that spoil this simple picture, once one 
goes beyond very large scales, as will also be seen in this paper. 
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It is worth pursuing how far we can understand RSD for the simple reason that RSD offer a 
unique way to measure growth rate of structure formation [3], and also can provide tests of dark 
energy models and general relativity [4-6]. Generically, if one had a good understanding of the 
nonlinear effects, RSD would be the most powerful technique for these studies because redshift surveys 
provide 3-dimensional information, while other methods, such as weak lensing, provide 2-dimensional 
information (or slightly more if the so-called tomographic information is used [7, 8]). The most 
problematic part of RSD studies are the nonlinear effects, which have proved to be difficult to model, 
and which can extend to rather large scales, making their modeling essential for using the RSD as a 
tool. 

In recent years several studies have been performed investigating these effects [9-12]. Some of 
these studies included galaxies or halos, [13-16]. Some of these methods use analysis and modeling 
based on perturbation theory (PT) [17], but none attempt to rely entirely on PT to explain all of the 
effects. Instead, they rely on ansatzes with free parameters, so that if the ansatz are accurate one 
can model the effects accurately. Separately, there have been several approaches trying to improve 
perturbation methods and to increase their ranges of validity [18-27]. All of these approaches adopt 
a single stream approximation, which we know breaks down on small scales inside the virialized halos 
and which is particularly problematic for modeling of RSD. 

The goal of this paper is to present a systematic PT approach to all of the lowest order terms 
contributing to RSD. Our goal is to identify which can be modeled well with PT, which can be mod- 
eled with extended PT methods mentioned above, and which require phenomenological additions to 
account for the small scale physics which cannot be modeled with traditional PT that does not include 
velocity dispersion. This approach is enabled by the recently developed distribution function approach 
to RSD [28] , which decomposes RSD contributions into separate correlations between moments of dis- 
tribution function. As such it allows us to investigate individual contributions to RSD and develop 
different PT or other approximation schemes for these terms. Whether this is ultimately useful for 
modeling RSD remains to be seen: our primary goal is to develop better physical understanding of 
dominant contributions to RSD. 

The paper is organized as follows: we begin in Sec. 2 by presenting a more detailed derivation of 
the angular decomposition of redshift space power spectra than given in [28]. We then use in Sec. 3 the 
perturbative methods to model the lowest contributing terms in this expansion, augmented by simple 
phenomenological models and/or beyond the lowest order contributions to improve the model when 
necessary. Results are also compared to numerical simulation measurements presented in [29]. We 
summarize and conclude in Sec. 4. In Appendices A, B, C, D we show some details of the calculations 
and write explicit forms of the terms contributing to the power spectra. 

For this work, flat ACDM model is assumed Ct m = 0.279, 0\ = 0.721, Q, h /Q, m = 0.165, h = 0.701, 
n s = 0.96, <Tg = 0.807. The primordial density field is generated using the matter transfer function 
by CAMB. The positions and velocities of all the dark matter particles are given at the redshifts 
z = 0, 0.509, 0.989, and 2.070, which are for simplicity quoted as z =0, 0.5, 1, and 2. 

2 Redshift-space distortions form the distribution function 
2.1 Generation of velocity moments 

Evolution of collisionless particles is described by the Vlasov equation [30] 

d f = ^- + — -V x f- amV<j> ■ V p f = 0, (2.1) 
cIt ot am 

where the gravitational potential (f> is given by 

V%<i> = ^Ga 2 p5 = ^n 2 n m S. (2.2) 

Here /(x, p, r) is the particle distribution function at a phase space point (x, p), where x is the 
comoving position, p is the corresponding canonical particle momentum defined by -jE = —am\7<p. 
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t = J dt/a is the conformal time, m is the particle mass, and % = dlna/dr — Ha is the conformal 
expansion rate, where H is the Hubble parameter. 

Note that in this paper we will use the canonical momentum p rather then comoving q = p/a 
defined in [28]. The reason is that the comoving momenta q is not the canonical coordinate to 
comoving position x, and this would lead to additional terms in the Vlasov equation (because of 
coordinate transformations), i.e. taking corresponding q- moments of usual form of Vlasov equation 
2.1 would not give the standard form of continuity equation, Euler equation, and higher moment 
equations. This is not a inconvenience when the symmetries are to be considered, but in order to 
avoid this we will use the canonical momenta p. 

In the following we will drop explicitly writing the time dependence, i.e we will write /(x, p). 
The density field in real space is obtained by integrating the distribution function over the momentum 
space 

p(x) = ma- 3 J d 3 pf(x,p), (2.3) 

and mean (bulk) velocity v and higher moment fields can be similarly obtained by multiplying the 
distribution function by corresponding number of particle momentum p = amu (u is here a particle 
peculiar velocity) and then integrating over it. The mean velocity field of a particles is then given by 

f d 3 pJ2-/(x,p) 

v x = J , ;™ M 7 , 2.4 

and the velocity dispersion tensor is 

^(x)^ /d r 3 ^ /(x ; P) -,V, (2.5) 
J d 6 pf[x, p) 

i.e. c^x) = (dv % dv^ with dv l being the deviation of a particle's velocity from the local mean veloc- 
ity, and the average is taken over all particles at position x. Note the difference between the particle 
velocity u and mean velocity v. The first one is the velocity of a single particle that corresponds to 
the canonical momentum p, which is one coordinate in the phase space. On the other hand v is a 
field defined at every coordinate x and is averaged over all the phase space. In the similar way higher 
order moments can also be considered. 

Taking a arbitrary constant unit vector h, we can construct a following object 

T h V)-— /^/(x, P )(^) L , (2.6) 
P J \ma J 

i.e. velocity moments projected on the direction of vector h, and where p is the mean mass density. 
If we introduce approximations in which we neglect velocity dispersion and anisotropic stress, i.e. 
we neglect all the contributions from this second rank stress tensor, and similar higher rank tensors 
(er lJ = 0, . . .) it can be shown (App. A) that 2.6 is reduced to 

T h i (x) = (l + <S(x))(h.v(x)) i , (2.7) 

where 8 is a usual overdensity field (5 = p/p — 1). 

In this paper we omit the following Fourier transform (J 7 ) conventions 

/» = T [/(x)] (k) - J d 3 x exp(*k • x)/(x), 

/d 3 k 
exp(-k-x)/(k). (2.8) 
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2.2 Redshift-space distortions 

In redshift space the position of a particle is distorted by its peculiar velocity, thus the comoving 
redshift-space coordinate for this particle is given by 

s = x + f^, (2.9) 

where r is the unit vector pointing along the observer's line of sight, uu is radial comoving velocity, 
amity = p\\ = p f. The mass density in redshift space is then given by 

Ps (s)=ma- 3 J d 3 pd 3 x f(x,p)6 D (s-x-f^) = ma' 3 J d 3 p f (s - r^, p) . (2.10) 

By Fourier transforming equation 2.10, we get 

p s (k) = ma- 3 J d 3 x d 3 P / (x, p) e (*-*+**ii«n/«) 

■J d 3 x e ikx f d 3 P f(x,p)e ik ^ u ^ n , (2.11) 



ma 3 



were k is the wavevector in redshift space, corresponding to redshift-space coordinate s. 
Expanding the second integral in equation 2.11 as a Taylor series in k\\u\\/mH, 



ma 3 



f d 3 P /(x,pK fc nv« = ma -3 r d3 g /(x,p)5^i(iA,| U ||/w) j 



(2.12) 



where in the last part we have used equation 2.6 setting the vector h to be the unit vector pointing 
along the observer's line of sight h = f. Using that in equation 2.7 we have, in the case with no 
velocity dispersion or other second or higher rank tensors (which we will not generally assume) 

T|f(x) = (l + 5(x)>f(x). (2.13) 

The Fourier component of the density fluctuation in redshift space is 



(2.14) 



L=0 

were Tjf(k) is the Fourier transform of the Tjf'(x). For L=0 we drop the unmeasurable k = mode, 
and we are left with the density fluctuation Tj?(k) = <5(k). 

2.3 Angular decomposition of the moments of distribution function 

In order to make the context of this paper more clear we repeat angular decomposition of the moments 
of distribution function from [28], providing more detailed derivation. The object (x) introduced 
in equation 2.6 can be obtained as taking all components of moments of distribution function in h 
direction, which are the rank L tensors, 

—3 r 

Tti 2 ,... iL = J d 3 pf(K,p)u llUl2 ...u lL . (2.15) 

The real-space density field corresponds to L = 0, i.e. zeroth moment 2.3, the L = 1 moment 
corresponds to the momentum density 2.4, L = 2 gives the stress energy density tensor 2.5 etc. These 
objects are symmetric under exchange of any two indices and have (L + 1)(L + 2)/2 independent 
components. They can be decomposed into helicity eigenstates under rotation around k. 
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The full detailed derivation of this decomposition is done in (App. B) and here we give the final 
result taking h = f; 

m—l 

T lf( k ) = E E nf2f' m (fc)*WM), (2-16) 

(l=L,L-2,...) m=-l 

where coefficients nf are defined in equation B.7, and spherical tensors T^' m in equation B.5, and 
evaluated in frame where z || k, so it does not contain any angular dependence. 

2.4 Redshift power spectrum 

In our analysis we will adopt a plane-parallel approximation, were only the angle between the line of 
sight and the Fourier mode needs to be specified. The redshift-space power spectrum is defined as 
(S s (k)S*(k')) = (2vr) 3 P ss (k)(5 £l (k-k'). Equation 2.14 gives, 

where (2%fP w (k)S D {k - k') = (l]f (kJTJ^'fV)). Note that P LL >{k) = P L 'L(k)* so that the total 
result is real valued (what is explicitly shown in PT approach in App. C). Thus only the terms 
Pi£/(k) with L < L' need to be considered, each of which comes with a factor of 2 if L ^ L' and 1 if 
L = L' . If we introduce /i = feii/fe = cos#, we can write, 

^> = E(^(ff ^>+^£ E ^(f <«*> 

Zj — .£/ — L r L 

Next we insert the helicity decomposition of equation 2.16 and consider the implications of rotational 
symmetry on the power spectrum. Each term Plv (k) contains products of multipole moments 



rpL,' 



Xk)Y lm (9,cj>) (if ' m '(k)ty, ro ,(M))* « e l(m -" 1 '^. (2.19) 



Upon averaging over the azimuthal angle 4> of Fourier modes all the terms with m ^ ml vanish. 
Another way to state this is that upon rotation by angle ip the correlator picks up a term e^ m ~ m w 
and in order for the power spectrum to be rotationally invariant we require m = m! '. Putting all these 
together we find 

l 

p w (k)= Yl E E p(f'' m ik)ir(i>)irQ*), (2.20) 

(l=L,L-2,...) (l'=L',L'-2,...) m=0 

where P™^ — cos 9), are the associated Legendre polynomials, which determine the 8 angular de- 
pendence of the spherical harmonics. We absorb all of the terms that depend on I and m and various 
constants into the definition of power spectra P^{^ ' m (k). Note once again that these spectra depend 
only on amplitude of k. We have 



Pi 



f' m {k) cx (^T l L - m (k) I (r^' m '{k)X\ . (2.21) 



We also replaced the two helicity states ±m by a single one with m > 0, since their 9 angular 
dependencies are the same, and we absorbed the appropriate factors into the definition of P^ ' m (k). 
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2.5 Perturbation theory approach 

The parameter of the expansion in equation 2.18 can roughly be defined as k/iv/'H, where v is related 
to a typical gravitational velocity of the system. This velocity should be of order of a few hundred 
km/s, but note that higher and higher powers of these velocities enter the series. The expansion series 
is expected to be convergent if expansion parameter is less then unity. 

The main goal of these paper is to use perturbation theory to compute and assess contributing 
terms in expansion formula 2.18 in next to leading order (one loop). There is a close, but not one to 
one, relation between the PT expansion and expansion in distortions function moments. Assuming 
that S and kv/Ji make the same order of expansion in one loop (except P04 where leading term is two 
loop quantity) regime we obtain 



P ss (k) = P 00 (k) 




P 22 (k) + 2Re 



P 03 (k) 



—ikfi 
H 

kfj, 

u 



P)i(k) 



Pi 2 (k) 



P) 4 (k) 



(2.22) 



Neglecting all the velocity dispersion and anisotropic stress contributions we can use simplified form 
of Tjf (equation 2.13). After preforming the Fourier transformation we obtain 



T | f(k)=7-[(l + 5(x))^(x)] (k). 

In one loop PT regime only first three momenta are needed, so we can write 

lfl(k) = «n(k) + («||oa) (k), 

I|f (k) = o (k) + (v { \ o Uy o S) (k), 



(2.23) 



2]?(k) 



where we have used the following convention for convolution 

d 3 q 



(fog)(k) = J ^/(q) 5 (k-q). 



(2.24) 



(2.25) 



From the approximations we have adopted it also follows that curls of velocity field can be neglected, 
i.e. V x v(x) = 0. Thus velocity field can be fully described by its divergence field #(x) = V • v(x). 

So it follows uy(k) = i^6(k). 

At this point it is useful to observe that if working in linear perturbation regime well known 
Kaiser result [1] can be easily obtained directly from equation 2.17. It follows 



P ss (k) = P 00 (k) + 2Re 



and after using the facts that P^ ' 1] (k) = -fHP L {k) and that Pge''{k) = (/"H) 2 P L (k) we obtain 



-ikfj, 



Poi(k) 



kjxY 

u) 

(1,1), 



Pn(k), 



(2.26) 



P ss (k) = (l + / M 2 ) 2 P 5(5 (k), 



(2.27) 



hence, the Kaiser formula. 
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3 Perturbation theory results and comparison to the N-body simulations 



All of the N-body results used here have been presented in [29] . Briefly, for all of the power spectra of 
the derivative expansion one needs mass-weighted velocity moments, which can be straightforwardly 
measured from simulations. In [29] a series of TV-body simulations of the ACDM cosmology seeded with 
Gaussian initial conditions has been used [31]. We employ 1024 3 particles of mass 3.0 x lO n /i _1 M0 in 
a cubic box of side WOOh^ 1 Mpc. We use 12 independent realizations in order to reduce the statistical 
scatters. For the details of the simulations measurements we refer to the [31] and here we shortly 
repeat the basics. 

3.1 Poo(k): the isotropic term 

At the lowest order in kv/T-L expansion we have auto correlation of density field TjP(k) = <5(k). Power 
spectrum, Po (k)S D (k — k') = (5(k)<5*(k')), is well known and has been intensively studied, e.g. 



oo I 

[17, 32]. This first term does not have any \i dependence since it is independent of red shift space 
distortions, it dominants for small values of /i and in the limit fi = the transverse power spectrum 
becomes overdensity power spectrum P 00 (k). On scales smaller than k~ x ~ IQMpc/h, nonlinear 
corrections increase the power over the linear. 

Familiar one loop PT result for overdensity power spectrum is [17] 



Poo(k,r) = P 44 (k,i-) 



D\t)P& 1 \\.) + D\t) [pjf } (k) 



2P 



(1,3) 



Art 



(k) 



(3.1) 



where we have restored time dependence, with D(t) being linear cosmological growth factor. Pgg' 1 ^ (k) 
is the linear power spectrum Pf, (k) , and one loop contributions are 



P 



(2,2) 



Art 



(fc) 



(2tt) 3 
2/oo(fc), 



P L (q)P L (\k~q\) k (s) (q,k-q) 12 



P^ 3 \k) = 3P L (k) J ^P L {q)F^ (k,q,-q) 
= 3k 2 P L (k)J 00 (k). 



(3.2) 



Explicit form of all integrals of the I mn (k) and J mn (fc) type can be found in App D. In figure 1 one loop 
PT results for power spectrum have been presented, along with some of the other approaches, such as 
the closure theory approach [24] obtained from the Copter code [32] and the semi- fitting method [33] , 
based on power spectrum obtained from Zel'dovich approximation. Note that if one wants to impose 
consistency in expansion 2.22 and PT approach, only one loop regime PT result should be considered. 
All the power spectra on the figures are divided by the linear power spectrum fitting formula from 
[34] without BAO wiggles. We see that none of the methods give perfect agreement across all range of 
scales. SPT (one loop PT) actually gives the best results for k < 0.05h/Mpc, but predicts too much 
power at higher k. 

3.2 P 01 (k) 

The next term to consider correlates the overdensity field T^(k) = S(k) and radial component of 
momentum density P^(k). This is the dominant RSD term sensitive to velocities. As we can see 

from equation 2.16, momentum density can be decomposed into a scalar (m = 0) T t ' and two vector 
(to = ±1) components T-j 1 ' 1 . Only the scalar part correlates with the density Tq'°, which is a scalar 

field. Thus only non-vanishing contribution comes from PQ , 1 1,0 (fc) oc ^PQ : °(fc) (^T^°(k)j V what 

gives the simple angular dependence 



Pu(k) = p; 



o.i 



\k)P?(v)=nP^-°(k). 



(3.3) 
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Figure 1. Poo(^) power spectrum term is plotted at four redshifts z — 0.0, 0.5, 1.0 and 2.0. We show linear 
result (black, dotted), one loop PT (blue, solid), two loop closure (green, dashed), corrected Zel'dovich (red, 
long-dashed) of [33] , simple Zel'dovich (magenta, dot-dashed) and simulation measurements (black dots). The 
error bars show the variance among realizations in simulations. The power spectrum is divided by no-wiggle 
fitting formula from [34], to reduce the dynamic range. 



On the other hand, correlating directly (<5(k)T]|' 1 (k', r) ), from the equation 2.24 one gets power 
spectra 



p 01 (k) = -i^p se (k)-iA 01 {k), 



(3-4) 



where the first term is also well studied correlation function of overdensity field and divergence of 
velocity field 



(2^) 3 A 01 (k)^(k-k')= J ^3^(5(k)r(q)5*(k'-q)>. 



(3.5) 



For the first term, correlation function of overdensity and divergence of velocity field, one loop 
PT gives 



Pse(k,r) = D 2 (r)P^>(k) + D\t) P^>(k) + 2P&*>(k) 



(3.6) 



where Pgg (k) = — j'H.Phik} is the contribution in the linear regime, and one loop contribution is 

P& 2 \k) = -2m J ^P L (q)P L (\k-q\)F^ (q.k-q)^ (q,k-q) 
= -2fHI i(k), 



P$ 3) (k) = -3fUP L (k) 



^ P dl)\ [f [ , s) (k,q,-q) + G< s) (k,q,-q) 



-3fHk 2 P L (k)J 01 (k) 
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_j_ p ee 3) ( k ) 

{my 



Plr\k) + tM^). (3.7) 



Here we have introduced logarithmic growth rate / = /(r) = din D/dlna. 

For the second term in equation 3.3, we expand all the fields to the second order, i.e., one loop 
in the correlation function. Schematically, this gives 

(60S) = (sW9W>) + (aW^W) + (8V-W>5<») , 

or in therms of power spectrum 

A 01 (k,r) = D\t) (4i U) (k) + 4i 12) (k) + 4i ia) (k)) ■ (3-8) 
Again, using one loop PT we obtain the contributions from each of the terms 



4f(k)=-2/^P £ (k) / ^ g_M^(-q,k)Ji( g ) 

= -2fHfikP L (k) 3J 10 (fc) + ^ (o* + |j 
< 2) (k)^2/HP i (k) | ^0 Mi?W (qjk )p i(?) 



/■«/'/<■/ l (AM K + |^] - f-'" 

where the ^ = | / ^ ^) is the one-dimensional velocity dispersion at linear order, and a 2 = 

J ^ ? 3 Px(g)' Note that all three terms give the same angular dependence, so Aoi ~ /i, and then 
follows that Pqi ~ /x, as was expected form the symmetry consideration on the beginning. Finally, 
collecting all the terms 3.4, 3.5, 3.8, 3.9 one loop PT prediction for the Pqi follows. Now the total 
contribution to the redshift power spectrum P ss from the Pqi term is 

P^(k,r) = 2=^P 01 (k,r) = 

= 2]{t)D 2 (t) [ i 2 (PL(k) + 2D 2 (t) [Ioi(k) + I w (k) + 3k 2 (j 01 (k) + J 1Q (k))P L (k)]) . (3.10) 

In this form result is naturally separated in linear and one loop contribution part. Note that linear 
part here is the second term of Kaiser formula. 

Alternatively, the scalar mode of momentum can be obtained from the divergence of momentum 
and related to S using the continuity equation 6 — ikp s = 0, which is in terms of quantities defined 
previously 

T °<° - ikT^° = 0. (3.11) 

Note that the vector part of momentum field does not contribute, since it vanishes upon taking the 
divergence (i.e., vector components are orthogonal to k and the dot product is zero). 
It follows 



i%i(k J T) = i-P j4 (A,r)=t— — , (3.12) 
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and the total contribution to P ss (fc,r) is 



P s f(k,T)= M 2 ^- 



dP Q o(k,T) 

dr 



2 «■« oo(fcj a) 



,dRu 



dint 



(3.13) 



This result, first obtained in [28], is exact for dark matter, valid also in the nonlinear regime. It shows 
that this term can be obtained directly from the redshift evolution of the dark matter power spectrum 
Poo(fc), so if we have an accurate PT model for P 00 then we should also have the same for P i- On 
large scales it agrees with the linear theory predictions. If we write Poo(fc, t) = Z)(r) 2 Pi(fc, r), we find 
Kaiser part Pfif = 2//i 2 P/j n (fc). On smaller scales we expect the term to deviate from the linear one, 
just as for P 00 (k). Using one loop PT we simply need to calculate the derivatives of growth factor 
D(t) = f(T)HD(r), and from equation 3.1 we get 



P^(k,r) = f(r)HD\r) [p^\k) +2D 2 (r) (p^ 2) (k) + 2P^ 3) (*)) 



Finely, plugging that in equation 3.13 we get 



P t(k) = 2/(r) J D 2 (r) M 2 



P L (k) + TD 2 (t) (loo (k) + ik 2 J m {k)P L (k)) 



(3.14) 



(3.15) 



After some integral transformations and calculations it can be shown that this result is equivalent to 
the on in equation 3.10. Obtained PJ-f results are presented in Figure 2. We show the one loop PT 
results, along with semi-fitting method [33] based on power spectrum in Zel'dovich approximation, 
and simulation measurements. The power spectra are now divided by second term in Kaiser formula 
where no-wiggle linear power spectrum has been used. 

3.3 P n (k) 

The next term we are to consider is the autocorrelation of momentum density ^(k) field. In this 



case scalar (m = 0) T i,u (fc) correlates with itself, and the vector (m = ±1) components T^' ±1 (k) also 
correlate with itself, so both components of momentum contribute, 

Pu(fc) = P 1 1 ; 1 1 '°(fc) [PfOu^ + P^fc) [PfV)] 2 - 

Contributions to redshift space power spectrum is then given with 



(3.16) 



(3.17) 



The scalar part is the autocorrelation of the of the momentum that contributes to the continuity 
equation 3.11. In linear PT only the scalar contribution in non-zero and P 1 1 ' 1 1 '°(fc) = / 2 P^(fc), which 
is the last term in Kaiser formula. There is another contribution to both /x 2 and /i 4 terms from the 
vector part of momentum correlator P 1 1 ' 1 1,1 (/c) oc (^T^' 1 [k)\ 2 \ , which comes in at the second order in 
power spectrum, and can be computed using one loop PT. This vector part is often called the vorticity 
part of the momentum, because vorticity of momentum does not vanish, even if vorticity of velocity 
vanishes for a single streamed fluid [35]. From equation 3.17 can be seen that this term always adds 
power to term and subtracts it in /n 4 term, but is combined with a positive contribution from the 
scalar part in /x 4 term. 

Now using expressions 2.24 we can straightforwardly expand the correlator in density 5 and 
velocity divergence 9 fields. In terms of power spectra we have 



P n (k) = ^P M (k) 



2|p n (k) 



Cii(k), 



(3.18) 



where we have introduced: 



(27T) 3 P ee (k)^(k-k') = <0(k)0*(k')} 
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Figure 2. fc-dependence of Pq* term of redshift power spectrum is plotted at four redshifts z = 0.0, 0.5, 1.0 
and 2.0. This therm has simple fi 2 dependence in all nonlinear orders. Here we show linear Kaiser result 
(black, dotted), one loop PT (blue, solid), corrected Zel'dovich (red, dashed) model from [33], simple Zel'dovich 
(magenta, dot-dashed), and simulation measurements (black dots). The error bars show the variance among 
realizations in simulations. The power spectra are divided by second, no-wiggle, term of Kaiser formula to 
reduce the dynamic range. 



(27T) 3 i? 11 (k)^(k-k') 

(2n) 3 C 11 (k)d D (k-k l ) 



^<%)<Kk-q)0*(k')>, 



(3.19) 



Using one loop PT to evaluate these power spectra. First term gives familiar velocity divergence 
autocorrelation 



P ee (Kr) = D(rfP^\k) + D\t) (p£ 2 \k) + 2P^(k)) , 



(3.20) 



where Pgg'^ik) is the linear power spectrum (fH) 2 Pi,{k) and rest is one loop contribution to velocity 
divergence power spectrum Pgg(k), 

P d i2 e 2) (k) = 2(fH) 2 [ ^Lp L (q)P L (|k-q|)[G«(q,k-q) 



/ (2tt) 3 
= 2(fH) 2 I 11 (k) 

P^' 3 \k)^3(fH) 2 P L (k) J ^Pl(i)G { 3 S) (k,q,-q) 

= 3(fH) 2 k 2 P L (k)J n (k) (3.21) 
Second therm can be expanded in the fields to the second order; schematically we have 

(689) = (ewswow) + (eWWv) + (e^s^e^ 
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This gives in terms of the power spectrum Pn(k) 

B u (k,r) = D 4 (r) (sg U) (k) + fl£ 12) (k) + b{\ 12) (k) 
where contributing terms are 



(3.22) 



B^\k) = 2(fU) 2 P L (V / ^ £^G«(k,- q )P i(q ; 



2(fH) 2 ^kP L {k) 



3Jio(fc) + 
d 3 c 



1 



3/fc 2 



< 1) (k) = 2( / H) 2 F L (k) / ^ ^i^-q^q 



-(fnYfxkP L (k) ( ^ 



3/c 5 



B[\ 12 \k) = 2(fH) 2 J ^ ^G 2 s) (k-q,q)P L (k-q)P L (q) 



= 2(fU) 2 ^I 22 (k). 
Similarly, for the last term in equation 3.18, we have 



(3.23) 



c[ 1 1 ni \k) = (fny 



d 3 q q\\ 



(2n)° q 



3 „2 



(k-q)| 
(k-q) : 



P L (k-q)Pi(q) 



(fH) 2 k- 2 (l 31 (k) + v 2 I 13 (k)) 



(3.24) 



Combining all that, we can write the contribution to redshift space power spectrum P ss from 
Pn term 



AT(k) 



Pn(k) = f 2 {r)D 2 { T )^ 2 (/. 2 P L (fc) +P 2 (r)/ 3 i(£ ; )) 



+ / 2 (t)D 4 (t)m 4 [2Ju(fe) + 4/ 22 (fc) + J 13 (fe) + 6fe 2 (J u (fc) + 2J 10 (fe))P L (/c)] 



(3.25) 



As can be seen we obtained /i 2 and /i 4 angular dependence from this term, as was argued from 
symmetry consideration in [28] . Vector contribution can be identified as the part multiplying /i 2 [28] . 

On the other hand, we could have started directly from equation 3.16. If we chose to work in 
the frame where z || k one can write the decomposition 2.16 of momentum density Tjf (k) = Py(k) = 
f • p(k) = p s cos 9+p v sin 9 cos <fr, where we have chosen, without loss of generality, for f to be in x — z 
plane, and p s and p v represent scalar and vector part of decomposition, respectively. After averaging 
over cf> angle, this enables us to write P^\' — Pp*,p B an d Pi'i' 1 = Ppv,pv- Just as before, scalar part 
can be determined directly from continuity equation 3.11. We can again use one loop PT to evaluate 
scalar and vector contributions 

Pp.*. = *- 2 Ps,t = (W) 2 D(r) 2 k- 2 (pi, M) (k) + D 2 (r) [4P] 2 ' 2) (k) + 6P^ 3) (k) 

= (fH) 2 D(T) 2 k- 2 (P L (k) + D(t) 2 [8I 00 (k) + 18k 2 J 00 (k)P L (k)]), 



P, 



d 3 q 



k x q 



2k q 



(2tt) 3 q 4 



(k - q)* 



P L (k)P L (\k-q\) 



= (fH) 2 D A (r)k- 2 I 31 (k). 
Thus, the contribution to the total red shift power spectrum from Pn term is 



Ptt=f 2 (r)D 2 (r)n 2 



fi 2 P L (k) + n 2 D(T) 2 (8I 00 (k) + 18J 00 (k)P L (k)) + (1 - A1 2 )P(r) 2 / 31 (fc) 



(3.26) 



(3.27) 
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Figure 3. fc-dependence of the scalar part of P"" term. Power spectrum is plotted at four redshifts z = 
0.0, 0.5, 1.0 and 2.0. This term has a simple fi 4 dependence. Here we show linear Kaiser (black, dotted) and 
one loop PT (blue, solid) result, and compare it to simulation measurements (black dots). The error bars 
show the variance among realizations in simulations. The power spectra are divided by the no- wiggle linear 
term. 

Again, after some coordinate transformations and algebra it can be shown that this result is equivalent 
to the one we obtained earlier in equation 3.25. 

In order to improve our prediction for the vector part we can take into consideration the most 
relevant higher order loop contributions. Starting from definition of C\\ (equation 3.19), which gives 
raise to the vector part of Pu, and generalizing our one loop prediction in equation 3.24 we can write 

Cll(k)= / ^(q)^(k-q) + / ^J^^ P MPse(k- q ). (3.28) 

In low k limit this gives back the previous result from equation 3.24, and in high k limit the first term 
dominates what is giving (fHD) 2 P^(k)a1 for the vector part. 

In figure 3 we show scalar part of P-f-f which comes from scalar contributions. It has a simple 
H A angular dependence, and corresponds to the third Kaiser term. We divide the plots by this Kaiser 
limit, using the no- wiggle linear power spectrum. One loop PT results are compared to the simulation 
measurements. We see that PT is quite successful in reproducing the nonlinear evolution of this term. 

In figure 4 we show the vector ^ 2 (1 — /z 2 ) part. We see that one loop PT is successful in reproduc- 
ing the simulations for k < 0.2h/Mpc (the disagreement for k < 0.03h/Mpc is likely numerical), and 
adding two loop corrections increases these rage to larger k. We also see that this vector contribution 
is considerably smaller than the scalar part for /i = 1 for most of the fc-range shown here, becoming 
comparable only at k ~ Q.hh/Mpc. However, because this vector term scales as /i 2 while the linear 
scalar term scales as /z 4 , the vector terms always dominates for sufficiently small fi. So for p, = 0.1 
the nonlinear vector part exceeds linear scalar part already at k ~ 0.05h/Mpc. 
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Figure 4. fc-dependence of scalar and vector part of P(f term of the redshift power spectrum is plotted at 
four redshifts z = 0.0, 0.5, 1.0 and 2.0, assuming fi = 1. Scalar part has simple fi 4 angular dependence while 
the vector part has fJ?(l — fJ?) angular dependence at all (nonlinear) orders. We show linear/Kaiser result 
(black, dotted), one loop PT result for scalar part (blue, solid), one loop PT result for vector part (lighter 
red, dashed), relevant part of two loop PT for vector part (red, solid) and simulations for scalar (blue points) 
and vector (red points) part. We also show scalar contributions of Cii term at one (lighter green, dashed) 
and two (green, solid) loop order. 



3.4 P 02 (k) 

At orders higher than Pn(k) there are no linear contributions, hence these terms are usually not of 
interest for extracting the cosmological information. However, these terms are known to be important 
on surprisingly large scales. These terms have usually been modeled phenomenologically in terms of 
adopting a simple functional form for k and /i dependence and are related to the so called Fingers- 
of-God (FoG) effect. We begin with the dominant P02 term, which, as we will see, is the last term to 
contribute to /i 2 dependence. 

We correlate the scalar density filed Tq'° = S with the tensor field T 2 ' m . Since scalars only 
correlate with scalars, there are only two different terms that contribute [28], 

P 02 (k) = P ° 2 -°(fc) [P»] 2 + P$°(k)P°(»)P°(n). (3.29) 
In terms of the contribution to the redshift space power spectrum this gives 



r Q2 



ir p^\k) + \p^(k)^-l) 



(3.30) 



The first term is the correlation between the isotropic part of the mass weighted square of velocity, 
i.e. the energy density T 2 ' = (1 + 5)v 2 , and the density field T ' =5, and the second term comes 
from the scalar part of the anisotropic stress T 2 1 correlated with the density Tq'° = 5. 

Before using PT to model these terms let us consider what we can expect from physical grounds. 
As argued in [28], in systems with a large rms velocity, the first, isotropic part Pq'q' should scale 
as P o(fc)<7 2 , where er 2 has units of velocity square and includes the small scale velocity dispersion 



- 14 - 



generated inside nonlinear halos. Some of this contribution cannot be modeled with simple fluid based 
PT, since not all of velocity dispersion is captured in this approach. As a result, we should not even 
hope that PT can be reliable for this term: we will need to add an extra contribution to account for 
the small scale velocity dispersion. 

Expanding the fields we can write the contributing terms as following 

P 02 (k,r) = -D 4 ( T )(A 02 (k) + B 02 (k)) ) (3.31) 

where we have 

(2 7 r)^ 02 (k)^(k-k')= / ^ ||^|(W(qr(k'-q]}, 

(2^) 3 P 02 (k)^(k-k') = J ^ ^(S(k)9*(q)e*( q ')5*Ck'- q - q ')}. (3.32) 

Using the one loop PT to evaluate this terms we expand these terms in the following way 

A 02 (k) = A& n \k) + A& 21 \k) + A& 12 \k), 

B 02 (k) = Bg ni \k), (3.33) 
which after some computation give 

= {fUfk- 2 (l 2(k)+n 2 ho(k)), 

4 1 2 21) (k) = 4 1 12) (k) = 2 ( /H) 2 p L (fc) J ^|^4 s) (q,-k)P,(q) 

= (mfp L {k) (j 02 (fe) + /i 2 j 20 (fc)) , 



Bg 111 \k) = -(fH) 2 P L (k) J ^3 ^PM) 

= -(fH) 2 P L (k)a 2 v . (3.34) 

Putting together all of the above we obtain for the Po 2 contribution to the total redshift power 
spectrum 



P s |(k,r)=-( ^) P 02 (k) 



/ 2 (r)D 4 (r)/x 2 



It 



, (k) + k 2 (2 J 02 (k) - a 2 ) P L (k) + /i 2 (/ 20 (fc) + 2k 2 J 20 (k)P L (fc)) 



(3.35) 



As we mentioned above, we have the contribution of form —/j 1 2 k 2 (j 2 Pi J (k), which suppresses the linear 
power spectrum with a k 2 like effect, increasing towards higher k. This is the lowest order FoG term, 
which we see contributes as (fc/1) 2 and so effects the fi 2 term of total P ss . Because small scale velocity 
dispersion effects cannot be modeled by PT, which is restricted to the weakly non-linear regime, we will 
consider a model where we add to the PT predicted value for velocity dispersion a 2 the contributions 
coming from small scales. In the equations above we can then replace a 2 — > a 2 + <7q 2 / (fHD) 2 , where 
<Tq 2 is the small scale addition to the velocity dispersion, and which we treat here as a free parameter. 
List of values used here for these parameters (depending on redshift), is given in the table 1, in 
section 3.8. In these section we also consider the explanation of these values using the halo model, 
see for example [36]. In addition to small scale velocity dispersion model, we also include the most 
relevant higher order PT terms. If we consider higher order contributions to (5\Sv 2 ) term we see 
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that it has subsets of diagrams where (6\ is not connected to any of the velocity fields, so we can 



write (S\Si 



(S\6) 



Formally, in one loop computation only the leading term of this subset 



contributes in equation 3.34. Collecting these we see that we can model Bq 2 term by replacement 



D\r)B m {k) = - (fHD 2 ) 2 a 2 v P L (k) -+ -(fHD) 2 (a 2 + a 2 2 /(fHD) 2 ) P 00 (k,r) 



(3.36) 



In order to discuss the results let us first rewrite equation 3.35 in form of isotropic and anisotropic 



part as for P 02 • We have Pq| = [i 2 \ P ( 

P{r)D\r) 



,SS,I 

02 



|(3/. 2 



i \ pss,A 
s -) Jr Q2 



where 



n S 2 J (k,r) = 



^or(k,r) 



3/ 02 (fc) + ho(k) + 2k 2 (3J 02 (k) + J 20 (k))P L (k) 

f 2 (T)D 2 (r)k 2 (a 2 + <r 2 02 /(fHD) 2 ) P 00 (k, r), 
2/ 2 (r)P 4 (r) \ + ^j^p^ 



(3.37) 



In figure 5 we show isotropic and anisotropic part of the Po 2 contribution to the total redshift power 
spectrum. All power spectrum contributions are divided by the (fDk) 2 a 2 P™(k), where we again 
used the no-wiggle power spectrum. We can see that the contribution to /j? is always negative, while 
the corresponding vector term from P n always adds power and the two partially cancel out [28]. As we 
see the scalar anisotropic stress-density correlator Pq 2 '° contributes to the fi 2 angular term, as well as 
to the /i 4 term. The anisotropic term is reasonably well modeled with PT and has smaller magnitude 
than the isotropic term, as expected, since the velocity dispersion in virialized objects is essentially 
isotropic. The isotropic term is poorly modeled with just PT: we need a significant contribution from 
the small scale velocity dispersion, which can be seen to essentially double the amplitude of this term 
at low k, and far more than that at high k. In figure 5 we can see that this term helps the model 
considerably, but of course it has one free parameter. 
We can also write this result in powers of /Lt, 

P$i\»*] = f{r)D\r)(l 02 (k) + 2fc 2 J 02 (fc)P L (fc)) - / 2 (r)P 2 (r)fc 2 (cr 2 v + <r 2 02 /(fHD) 2 ) P oa (k,r), 
= J&V] - f(r)D 2 (r)k 2 (a 2 + <x 2 2 /(/HP) 2 ) P o(fc,r), 



Po2\A = P&\A = f 2 (r)D\r)(l 20 (k) + 2k 2 J 20 (k)P L (k) 



(3.38) 



where we have implicitly defined Pq| by omitting the velocity dispersion part from P ( 



02 • 



3.5 



Pi2(k) 

As we can see the lowest order in [i with which correlators contribute is the L + L' or L + L' + 1. 
So contributions to /i 2 comes only from terms P01, Pn and Pq 2 . The next order in powers of /i 2 will 
be /i 4 terms. As we have seen P\\ and Pq 2 also have contributions to /i 4 , with P\\ having the linear 
order term which dominates on large scales. 

Here we correlate the momentum filed T \' m with the tensor field T 2 ' m . Because of rotational 
invariance we can correlate only scalar to scalar field and vector to vector field 



p 12 (k) = p 1 1 ; 2 '°(fc) [p8Qi)P?M] +P};i'°(k) [p?0i)P§M] +P 1 1 ; 2 2 ' 1 (fc) [Pfo)Pfoj\ 

In terms of the contribution to the redshift space power spectrum this gives 



P£(k) = -i 



kfj, 
H 



A,b ,U (*)A* + 2 ( fc )^(3M 2 - 1) + SP^'^Ml - A*") 



(3.39) 



(3.40) 



Using the one loop PT we get both /i 4 and /i 6 angular terms, but since there are 3 terms we 
cannot distinguish between them in equation 3.40. Using equation 2.24 and one loop PT we get 



P 12 (k) = -D*(T)p (k\\A 12 (k) + k 2 B 12 (k) + fc||Ci 2 (k)) , 



(3.41) 
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Figure 5. fc-dependence of isotropic and anisotropic part of Pq| term of redshift power spectrum is plotted 
at four redshifts z = 0.0, 0.5, 1.0 and 2.0. Isotropic part Pq^' 1 , computed in one loop PT (red, solid) 
is plotted, as well as using the model presented above (blue, dot-dashed). Isotropic part has simple fj 2 
angular dependence while the anisotropic part Pq2' A (g reen i dashed) has ^i 2 (3/^ 2 — l)/2 angular dependence. 
Simulation measurements (dots) for the corresponding terms are also presented. The power spectra are divided 
by fc 2 cr 2 P£ w without the wiggles. 



where the contributing terms are 
(2w) 3 A 12 (k)6 D (k-k') = J 



d 3 q qu (k' - q)i 



(2tt) 3 q 2 (k' - q) 



<0(k)0*(q)0*(k'-q)) 



d 3 qd?q' q\\ q\\ (k'-q')|| , , , , > 

W~ ^(k^F (#(k " q)fl {q)9 (k " q)) 



(27r) 3 C 12 (k)<5 D (k-k') = J 



d 3 q d 6 q' q\\ % 



<0(k)r(q)0*(q')^(k'-q-q')> 



(2tt) 6 g 2 q n 

The first of these terms we can be expanded further 

A 12 (k) = 4 2 2 n) (k)+4 m) (k)+4 U2) (k), 

and computing these terms gives; 



.4 



(211) 
12 



d 3 q q\\ (k - q)|| 



(k) ^- 2(/H)3 y (2.)3 g 2 (k-q) 2 

-{.fUfk- 2 (l 12 (k) + fi 2 I 21 (k)) 



G^ ) (q,k-q)P i (q)P L (k-q) 



42 ai) (k) = ^i 2 12) (k) = -2(fHfP L (k) J 



d 3 q q\\ (k-q)y^ (s) 



(2tt) 3 g 2 (k - q) ; 



Gr(q,-k)P L (q) 



-(/H) 3 Pi(k) (Jo 2 (fc) + /i 2 J 2 o(A;)) 



S ia (k) 



-2(/H) 3 



d 3 q 9M(k-q),| 



(2tt) 3 q 4 (k - q) 2 



P L (q)P L (k-q) 



(3.42) 



(3.43) 
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= (fH) 3 fik- 3 (l 03 (k) + ^I 30 (k)), 

C la (k) = (fH) 3 P L (k) J ^P L (q), 
= (fn) 3 a 2 v P L (k). 

All this gives us the contribution to total redshift space power spectrum 



(3.44) 



->\ ^) 7^>(k.r)=/lrri)(r)V' 



/i 2 (fc) - Jo3(*) + 2k z J m {k)P L {k) - k 2 a 2 P L (k) 



+ V 2 (hi(k) - ho(k) + 2fc 2 J 20 P L (fc)) 



(3.45) 



We can again add the small scale velocity dispersion in by hand, as was done and explained in case of 
Pq2- Considering the relevant higher order contributions we see that the isotropic part of (t^\Sv 2 
can be modeled by 



-iD 4 (T)|c 12 (fe) - -iD\r)^(fH) 3 a 2 v P L (k) —(fHD) 2 (<J 2 v +« 2 12 /(fHD) 2 ) P i(k,r), 

where we again treat small scale velocity dispersion a 2 2 as a free parameter with values for different 
rcdshifts given in table 1. These values are the same as for P02 case, and the reasons and explanation 
in term of halo model is given in section 3.8. 

As mentioned earlier, since we have only /j, 4 and /i 6 angular dependence we can not determine all 
three terms in equation 3.40 separately. Let us instead separate the angular dependences itself and 
collect the terms 

/':-;><': \{'- :r <> '':-; It 1 

= /(t) 3 £>(t) 4 [l 12 (k) - I 03 (k) + 2k 2 J Q2 (k)P L (k) 

\f{r) 2 D{r) 2 k 2 (a 2 v + a 2 12 /(fHD) 2 ) fft"(fc, r), 
= PS [|U 4 ] - \f{rfD{rfe (a 2 v + v 2 12 /(fHD) 2 ) P^(k,r) 



Pf| 



f(T) 3 D(T) 4 \l 21 (k) - I 30 (k) + 2k 2 J 20 P L (k) 



(3.46) 



where we again implicitly define Pf|, by omitting the velocity dispersion part. 

In figure 6 we show Pf| [/1 4 ] and Pf| [m 6 ] parts to the total redshift power spectrum. Power 
spectrum contributions are divided by the (fDk) 2 a 2 Pf w (k), where we again used no-wiggle power 
spectrum. 

3.6 P 22 (k) 

Next we consider correlator of tensor T 2 m field with itself. This term will give i/ 4 , /i 6 and /i 8 
contributions. One loop PT gives first order contributions to all of these angular terms. From the 
expansion of power spectrum 2.20 we can see that the constant contribution to P 22 is coming from 
the scalar term, Pq'q' and partially from Pq 2 '° and P^' - This will give /z 4 as the lowest order 
contribution to the total P 22 , and all of the other terms will come as n 6 , /z 8 . Let us now assess these 
contributions using one loop PT 



P 22 (k,T) = 2(fHD f J 



QJ (k-q)| 
Q 2 (k - q) 5 



P L (q)P L (k-q) 
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Figure 6. fc-dependence of fi 4 and fi 6 part of P£f term of redshift power spectrum is plotted at four redshifts, 
z = 0.0, 0.5, 1.0 and 2.0. /z 4 part Pol[/U 4 ] is computed in one loop PT (red, solid) regime, and using model 
presented above (blue, dot-dashed). We show /i 6 part Pol[/z 6 ] computed in one loop PT (green, dashed), 
and simulation measurements (dots) for the corresponding terms. All power spectra are divided by fe 2 <r^P|) w 
without the wiggles. 



= {.fHDf\k- A 



I 23 (k) + 2 f i 2 I 32 (k) + ^I 33 (k)), 



which gives rise to the total red shift power spectrum contribution 



^2 S |(k,r) = 



1 f k/i\ 4 



1 



H J 



P 2 2(k,T) = -/ 4 (r)£> 4 (r)/i 4 I 23 (k) + 2^1 32 {k) + ^I 33 (k) 
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(3.47) 



(3.48) 



These are the leading order contributions to the angular dependence of this term. Now let us also 
investigate the most important contributions from the higher orders. For that purpose let us write 
the full correlator in terms of the density and velocity fields 



In 



+ 2(v 



Svi 



(3.49) 



In equation 3.48 we have considered only the first of these three terms, but we should also include 
some of the most important contributions from the remaining terms. From two loop considerations 
first we improve the first term 3.47 by exchanging linear power spectrum Pl with one loop Pee- The 
most important contributions of the other two terms can be modeled as 



Sv \\ 
5vl 



(fHD) 2 a 2 v P 02 , 
(fHD)\o*) 2 P a 



P?2 ° P[){ 



These are of course not the only higher order term, but after a detailed analysis these terms turn out 
to be the most relevant and the rest can be neglected. We can again include the small scale velocity 
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Figure 7. fc-dependence of /i 4 , /i 6 and /i 8 parts of P|| term of redshift power spectrum is plotted at four 
redshifts, z = 0.0, 0.5, 1.0 and 2.0. fi 4 part P|| [t^ 4 ] is shown for one loop PT (red, solid, thin) regime, and for 
improved two loop PT model with small scale velocity dispersion (red, solid, thick), as well as for simulations 
(red dots). Similarly, Pf| [/x 6 ] part is shown using one loop PT (blue, dashed, thin) result, improved two loop 
PT model with small scale velocity dispersion (blue, dashed/dotted, thick), and simulations (blue, empty/full, 
diamonds). Dashed/full results present positive values of /i 6 dependence and dotted/empty negative values. 
/i 8 contribution is shown for one (thin, green, dot-dashed) and two loop PT (thick, green, dot-dashed), with 
the corresponding simulations (green stars). All power spectra are divided by fc 2 cr^P£ w without the wiggles. 



dispersion extending a 2 — > a 2 + a 2 2 j (fHD) 2 as we did previously for Pq 2 and P12. Combining all we 
obtain a model 

P 22 (fc,/i) =P m {k,n) - 2(fHD) 2 {a 2 v +a 2 22 /(fHD) 2 ) Pfofan) 

+ [fUDf (a 2 + o 2 2 /{fUD) 2 ) 2 Poo(fc) + (P22 o P 00 )(fc). (3.50) 

In high k limit last (convolution) term corresponds to 2(fHD) 4 a 4 Poo(k). In figure 7 we show the 
individual angular contributions for one loop PT calculus and for the improved model suggested above, 
and compare them to simulation measurements. We see that using the proposed model improves 
results in comparison to the one loop PT contributions, but still only qualitatively agrees with the 
simulations. One would find much better agreement if not imposing a 22 = o§ 2l i.e. with more free 
parameters. We mention that most of the correction to the /i 4 term comes from the isotropic modeling 



of the last ( Sv 



Sv 2 ) terms. Term 



also contributes to \i but less than the previous term. 

Corrections to the /1 6 come from the angular dependency of A§ 2 term and we see that it can explain 
the change of sign and scale growth trends. The additional terms do not affect the /z 8 term, which is 
well predicted (at least relative to /i 4 and with two loop PT model of first term. 

3.7 P 03 (k), P 13 (k) and P 04 (k) 

Our goal is to consider all terms at the ^ order. There are 3 left. First we consider terms P03 and 
P13. We correlate overdensity field or momentum field with rank three tensor field 2~j?(x). Angular 
decomposition for P03 is relatively simple since it has only scalar contributions, but P13 has scalar 
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and vector contributions. Using the angular expansion we get the following angular dependence 



p 03 (k) = Pg;l>\k)» + P^ik) 1 -^ (5^ 2 - 3) , 

P 13 (k) = Pt;?> {k)f + Ptf'Hk) (1 - fJ. 2 ) + Plf\k)\^ (5m 2 - 3) - Plf\k)l (5 M 4 - 6/i 2 + 1) 



(3.51) 



In one loop PT these terms are 



PM=lZ{fUfD\T)^P L {k) J 0^ P L (q) 
= i3(fH) 3 D i (T)vk- 1 P L (k)a 2 v , 

P 13 (k)=3(fHD(r))^ 2 k- 2 P L (k) J ^ (J) P L (q) 

= 3(fH) 4 D 4 (T)fi 2 k- 2 P L (k)a 2 v . (3.52) 

From angular decomposition of P03 we have scalar terms, Pq '- 3,0 and Pq 3' , contributing with angular 
dependence /i 4 and fj, 6 . One could evaluate these terms in PT, but at least two loop order is required 
for /1 6 , since in one loop order gives just p. 4 dependence. For Pi 3 we see that the lowest angular 
dependence comes from the vector contribution and not the scalar, although the scalar has lower 
perturbative order. Similar case was discussed for Pn, where the vector part, which is of one loop 
order, contributes to /j, 2 , while the leading linear order scalar part has /j, 4 dependence. 
One loop PT contribution to total P ss give 



PoU^r) = - P 03 (k) = -f A {r)D 4 {T)p 4 k 2 a 2 v P L {k), 

P-(k,r) = -1 P i3 (k) = ~f 4 (T)D 4 (T)n 6 k 2 a 2 v P L (k). (3.53) 

We can again include some higher order terms based on small scale velocity dispersion type arguments. 
For example, let us asses contributions to each of the terms above as if fully coming from 



Tf) = 3d 2 {T° 
T$) = 3(t 2 (t} 



T 1 " 



if 

and neglecting other two loop contribution. Taking this into account we can model terms above by 
replacing P 03 — > S(fHD) 2 a 2 P i and P13 — > 3(fHD) 2 a%Pii. In figure 8 we show results for both 
/i 4 part of Po 3 and for /i 6 part of P13. On the same plot we show one loop PT prediction for both 
terms (keeping in mind that in the overall contribution they differ relative to each other by the 
factor of — n 2 f(r). We compare model results presented above to simulations. The specific shape in 
simulations is explained by the proposed model, while it is not in one loop PT result. This effect 
arises from substitution of Pl with Pqi or Pn. We can again add the small scale velocity dispersion 
a 2 — > cr 2 + <TQ 3 /(fWD) 2 (or equivalently a 2 3 ) which is not included in PT analysis. In the model 
we suggest for Pi 3 we get, in addition to the /z 6 dependence, also \i 4 dependence, which comes from 
the vector part of P\\. In figure 9 we show that this can explain the trends and amplitude seen in 
simulations for this term. In the case of the leading contributions Po3[/i 4 ] and P 13 [/i 6 ] we found that 
lower value for small scale velocity dispersion cto 3 an d &13 is needed (table 1). This can be describer 
using the halo model and we return to that in section 3.8. Note that this value only affects the total 
amplitude, i.e. translates whole result up and down, but does not affect the shape. 

In a similar fashion we can estimate contribution of P04 term, which is the last term we need to 
consider at /i 4 order. Formally this term does not even contribute at the one loop order in PT, but we 
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can do two loop considerations as we did before. Considering the most relevant two loop contributions 
(from partially disconnected diagrams) this term can be modeled as 

(if I Ijf) = &{fHDfal(5 \vf) + 3{fHD)**$ (5 \S) + (S \S) («f \vf) c , 

where we used subscript c to label the connected part of the correlator. Here we again include the 
small scale velocity dispersion using a 2 — > crj; + <7q 4 / \fHD) 2 , just as in P 02 case. We can write the 
PqI contribution 



Pol [M 4 ] = - - 2 f{rfD{rfe {ol + alJ{fUDf) P™ [^} 

+ - A f{rfD{rYk i (a 2 v +a 2 J(fHD) 2 ) 2 P^(k,r) + ^P^(k,r) [ ^P 22 ( q , T ) 



Pot [/] = - \f{rfD(rfe {ol + *lJ{fUDf) j& [//] 



(2tt) e 



(3.54) 
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Figure 8. /i 4 dependence of p)| and /i 6 dependence of Pfl is plotted at four redshifts, z — 0.0, 0.5, 1.0 and 
2.0. One loop PT result is plotted (red, solid), as well as results of improved model discussed in the text for 
-Pol [m 4 ] (blue, dashed) and Pi|[/i 6 ] (green, dot-dashed). Results are compared to the simulation measurements; 
PoilMl ( blue > stars ) and P?£\M e ] (g reen > dots )- A11 the P lots are divided by no-wiggle -f{Dfi) 4 k 2 a 2 PE w for 
P S | and ffi 2 times this for the Pf| term. 



3.8 Halo model and small scale velocity dispersion 

In our analysis of correlators that contribute in expansion of the RSD power spectrum we find that 
some of them have terms proportional to velocity dispersion of dark matter particles. Particles moving 
in the gravitational potential can have large velocities even on the small scales, so they have significant 
contribution to the total velocity dispersion. Using PT we can evaluate velocity dispersion 

{fnD(T ^ = 3 J ( 27 r)3 ^ • (3 ' 55) 
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Figure 9. /i 4 dependence of Pf! and PqI term is plotted at four redshifts, z = 0.0, 0.5, 1.0 and 2.0. Simple 
modeled results for Pfl (blue, dashed), and PqI (green, solid) are shown and compared to corresponding 
simulation measurements; Pf| (blue, stars), and PqJ (g reen , dots). Both, simulations and model results for 
Pf| have negative values. All the plots are divided by no-wiggle (/_D/i) 4 fc 2 a^P£ w . 



Linear theory gives o~ v ~ 600fcm/s, but this does not properly take into account small scale con- 
tributions, which come from within virialized halos where PT cannot be used. Thus to take into 
account all nonlinear contributions we have to add this to our model in order to match the simulation 
predictions. In table 1 we show the values used for small scale velocity dispersion for modeled terms. 
These values are obtained by fitting our PT model using the free parameter for small scale dispersion 
(ace, 012; • • •) m order to match simulation predictions. We see that we can classify these into a few 
groups which have approximately the same value. 



Table 1. Small scale velocity dispersions as described in the paper (in km/s). 
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We can understand the fact that some terms have equal velocity dispersion to others, and some 
do not, using the halo model [36-40]. We can distinguish between three types of contributions to 
velocity dispersion. For terms P02, P12, -Pi 3 (vector) and P22 we find that the same value is needed. 
In these terms ujj always comes weighted by 1 + 6. In a halo model we divide all the mass into halos, 
such that the integral over the halo mass function times mass gives the mean density of the universe, 

JdM^M = p, (3.56) 

where dn/dM is halo mas function. Each halo has a bias b(M), which describes how strongly the 
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halo clusters relative to the mean, since 

1 f dn 

S(k) = -J dM — Mb{M)6(k). (3.57) 

Each halo also has a small scale 1-d velocity dispersion v 2 (M) oc M 2 / 3 , where the latter relation is 
only approximate and does not take into account effects such as halo profile dependence on the halo 
mass etc. 

We now decompose the terms into halos of different mass, accounting for small scale velocity 
dispersion v 2 (M), and accounting for biasing whenever this is multiplied by density S. For example 
for term P02 schematically we can write 

P02 = (S\(l + «5)«jf) ~ (*|*«f ) - (S\8) -_ J dM^Mb{M)v\{M) = P 00 < 2 , (3.58) 

i.e. we find that the velocity dispersion is weighted by bias. Note that we should have written 
the term 5 in halo model as well, but since the bias integrates to unity (equation 3.57) we do not 
have a contribution from the left hand side. Note also that we only include the small scale velocity 
dispersion effects here that come on top of the PT calculations above. Same quantity enters also in 
P\2, P\^{vector) and P 22 . 

For terms P03 and Pi 3 (scalar) we have a different contribution to small scale velocity dispersion 
because one of the velocity field in correlates with the density field and we can approximate 1 + 6 
with 1 at the lowest order. As a result v 2 is not density weighted. For example for P03 we have 
contributions from term 

P 03 = (5\(1 + S)vf) ~ (flvf) =3(*|«,|) -_ J dM^Mv 2 (M) = ZP Sv ^ 2 v2 . (3.59) 

Since there is no biasing and since b{M) > 1 at high mass halos which dominate the velocity dispersion 
these terms have a smaller value of velocity dispersion than we had in the first case. This is precisely 
what we find when fitting to the simulations. 

Finally, for the term P04 we find contribution 

P 04 = (<5|(1 + S)v\) ~ (s\8v\) = (5\S) -_ J dM^Mb(M)vl s (M) = P 00 a 2 bvi . (3.60) 

This term gives a value bigger then previous two because higher mass halos give a larger weight and 
they are more biased, which is also consistent with what we observe in simulations, and is presented 
in table 1. To convert v\\ into velocity dispersion we use the relation 

vl(M) = (235km/s) 2 , (3.61) 

see for example [41]. We use standard Sheth-Tormen model for halo mass function and halo bias [42]. 
We see that predictions from the halo model presented in 1 agree qualitatively but not quantitatively. 
This could be a consequence of the simplifying assumptions, such as ignoring the internal structure 
of the halo and its mass dependence. Note also that there are no errors in the analysis: it is possible 
that the sampling variance errors are large, specially for er 2 ^, which receives dominant contributions 
from the very high mass halos which may or may not be present in our simulations, depending on 
the realization. We do not go into a more detailed modeling here, but it is possible that with a more 
detailed model the agreement would improve. Even at this level the halo model gives an insight in 
hierarchy of the contributions er 2 2 < <j 2 v2 < o 2 ^ , and offers a qualitative picture why different terms 
in expansion need different values for velocity dispersion. 
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3.9 Putting it all together: /j, 2 - 5 terms, finger of god resummation and Legendre moments 

There are a finite number of velocity moment terms at each order of \i 2 ^ , in contrast to the Legendre 
multipoles expansion (monopole, quadrupole, hexadecupole etc), which receive contributions from all 
orders in moments of distribution function. We will thus investigate \i 2 i expansion, with the lowest 
3 orders containing cosmological information, while the rest can be treated as nuisance parameters 
to be marginalized over. Even in that case a good prior for these higher order angular terms would 
be very useful, although given the large number of terms that contribute to it it seems easier to be 
guided by the simulations rather than the PT. In this section we collect all the previous terms with 
fi 2 and /x 4 dependence. At \i 2 level the contributions come from Poi, -P02 and Pn terms, and for /x 4 
from Pn, P 2, P12, P22, P13, Pi 3 and P04 terms. In figures 10 and 11 we show xx 2 and /x 4 dependence 
of these terms divided by the corresponding no-wiggle Kaiser term. We show both the simplest PT 
model and the improved model that includes velocity dispersion effects. For modeling some of the 
terms we have been using the model for velocity dispersion a 2 — > a 2 + <t 2 j , where the added value a 2 ^ 
for term Py is given by the table 1. These model was optimized to fit corresponding terms primarily 
on large scales, where the dominant contributions comes from Pqi for /x 2 and Pn form /x 4 terms. To 
improve the model further for P i and scalar part of Pn, instead of PT predictions, we use exact 
values obtained from the simulations. We expect that ongoing activities in the modeling of nonlinear 
power spectrum will result in a successful model of these terms (note that P01 is given by the time 
derivative of the nonlinear power spectrum Poo). Although we have introduced a free parameters in 
our model note that P01 and Pn terms do not contain any free parameters, so we can use simulation 
results as well as any other method to predict these terms. 

The leading order in RSD is the /i 2 term. On large scales it is given by the Kaiser expression, 
but note that the deviations from the linear theory are of the order of 10% at z = already at 
k ~ 0.05h/Mpc. These nonlinear effects are dominated by the small scale velocity dispersion effects, 
which cannot be modeled by PT (a smaller effect, of the order of 2% at these scales, is caused by 
nonlinear effects in P01 which are modeled in PT) . This is a serious challenge for the RSD models and 
the ability to extract cosmological information from RSD: any additional free parameter that needs 
to be determined from the data will reduce the statistical power of the data set. Note however that 
we do not observe dark matter, but galaxies, so to address this concern in a proper way one will need 
to repeat this study with galaxies. We plan to pursue this in the near future. At higher redshifts 
these nonlinear effects are smaller: at z = 1 the 10% nonlinear suppression happens at k ~ O.lh/Mpc. 
In all cases the dominant nonlinear effect is to suppress the small scale power, as expected by the 
phenomenological models like [9], where a Gaussian smoothing is added to the extension of Kaiser 
formula. 

The fx terms show considerably more structure in the nonlinear effects: the overall power is 
initially suppressed relative to the linear term, stays flat for a while and then increases again (above 
k ~ O.lh/Mpc at z = 0). The effects are large: 20% suppression of power at k ~ 0.05h/Mpc for z — 
relative to linear. The model has some success in predicting some of these details, but is far from 
perfect and again it relies on the free parameters. The nonlinear effects are smaller at higher redshift, 
but remain significant. These /i 4 terms have an important contribution to RSD. For example, at 
higher redshift (where / ~ 1) they contribute about 30% to the quadrupole on large scales, with the 
dominant 70% contribution coming from /x 2 term. As for the /x 2 term, it remains to be seen how well 
we can model these terms such that we can extract the maximal information from the data, but the 
fact that the nonlinear effects are so large already on very large scales is a cause for concern. 

Our result for the RSD power spectrum P ss can be compactified in so called finger of good 
resummation. Following the ideas presented in [29] we can show explicitly that up to \i order our 
result in equation 2.22 can be written in the following way 

P£ G (k, fi) = exp [ - £V<g] (MQ + V 2 B(k) + n*C(k) + ...), (3.62) 

where we have defined 



A(k) = Poo, 



z=0.0 



z=0.5 




Figure 10. fi 2 dependence of P sa at four redshifts, z = 0.0, 0.5, 1.0 and 2.0. We show separately contributions 
of the PT (dashed lies), the improved velocity dispersion model (solid lines), and simulation measurements 
(points). The leading term is Pqi (red), to which we add Pq| (green), and to which we add P** to get the 
total (blue). Kaiser jj, 2 term (black, dotted) is also shown. All the lines are divided by no- wiggle fi 2 Kaiser 
term. 



B(k) = P ss [ t x 2 ]+k 2 a 2 oG P 00 , 

C(k) = F-[ M 4 ] + fcVf oG P s V] + \(k 2 <j 2 FoG ) 2 P m , (3.63) 

i.e. there is no change to fi 2 , and fi 4 terms, since the terms P ss [/j, 2 ] and P ss [/i 4 ] contain all the terms 
discussed previously and the FoG terms cancel by construction of B(k) and C(k). 

If we set the value of <J% oG = {fD) 2 a 2 this reduces to simple form where a 2 is now present only 
in the exponent in equation 3.62 and not in the brackets, i.e. 

A(k) = P 00 , 

B(k) = P^ + P^ 2 }+P^% 

C(k) = P>> [ M 4 ] + P S | [//] + P£ [//] + P^ [A , (3.64) 

where all the P ss and P ss terms here, as defined in previous sections, do not contain velocity dispersion 
contributions. This argument also generalizes to the case where we replace <Tp oG — (fD) 2 a 2 + cr 2 , 
where a 2 is the small scale velocity dispersion. This is the basic justification for using the FoG model. 

Unfortunately, the fact that a 2 cancels out is of limited use, since in practice the velocity disper- 
sion is not dominated by linear a 2 , but by small scale velocity dispersions, and as argued above, there 
is no single cr 2 , but instead there are several different velocity dispersions entering in the detailed 
RSD analysis at /x 2 and [i 4 order, cr bv 2, a v i and a^i. In fact, in our analysis we include these terms 
already so one can argue that it is the next term that we do not include that should enter in crp oG - At 
/j, 6 order there are again several velocity dispersions that can be defined and that have a wide range 
of values, so we cannot simply write down a value without explicitly evaluating all the terms at this 
order. It is however likely that their values will be of the same order as a v i and <J(,„4. In table 1 
we compare the root mean square average of these velocity dispersion values to the best fit value for 
(Tp Q, showing that indeed the value of <Jp oG is indeed related to these other values. 
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Figure 11. /1 4 dependence of P sa at four redshifts, z = 0.0, 0.5, 1.0 and 2.0. We show the PT (blue, dashed 
line) results, the improved velocity dispersion model presented in this paper (blue, solid line) , and simulation 
measurements (blue, stars). Also, the leading P"f term is shown in red. Kaiser /1 4 term (black, dotted) is also 
shown. All the lines are divided by no-wiggle /1 4 Kaiser term. 



It is customary to expand the redshift-space power spectrum in terms of Legendre multipole 
moments. The motivation for this is that when using the full angular information Legendre moments 
are uncorrelated on scales small relative to the survey size. Using ordinary Legendre polynomials 
Vi(n), we have 

p ss (k,ri= Yl pt'wnw, (3.65) 

(=0,2,4,-- 

where multipole moments, -P ; ss , are given by 

Pf s {k) = (2l + 1) [ P"{k,n)Pi{n)dn. (3.66) 
Jo 

where Vi{p) are the ordinary Legendre polynomials, Vo(fi) = L ^Mm) = (3^ 2 — l)/2 and V±{ij) = 
(35)U 4 — 30^i 2 + 3)/8. In the RSD analyses we are usually limited to modeling the monopole (I = 0) 
and quadrupole (I = 2) terms, although some information is also contained in hexadecapole term 

a =4). 

In figures 12 and 13 we show monopole and quadrupole power spectra predictions of improved 
velocity dispersion model presented in the paper, as well as one loop PT result. We also show 
resummed FoG result choosing for erpoG values given in the last line in table 1. We compare this 
to the reference multipole results obtained from full simulation redshift space power spectra. We 
also show simulation results where only terms up to /1 4 are considered. In case of monopole we see 
that these two simulation results agree on scales larger then k ~ (0.15 — 0.20)h/Mpc (depending on 
redshift) but then start to deviate one from an other. In the case of the quadrupole these deviations 
start to be more then 1% for k > 0.15h/Mpc. This trend is due to the /i 6 term which is weighted by 
1/7 for the monopole but 11/21 for the quadrupole, which is almost the same weight (4/7) as for /i 4 
term. At higher k higher [i terms (/i 6 , /i 8 , ...) start to be relevant and contribute significantly to the 
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Figure 12. Monopole moment Pq s is plotted at four redshifts, z — 0.0, 0.5, 1.0 and 2.0. To the first 
isotropic term (red) of P sa expansion we first add fj, 2 (green) term and then also fi 4 term (blue). We show PT 
(dashed lines) results, improved velocity dispersion model presented in this paper (solid lines), simulations 
up to (i 4 contributions (triangles, diamonds and stars), and reference simulation results (points). Resummed 
FoG model from equation 3.62 (black dot-dashed line) and linear Kaiser model (black dotted line) are also 
shown. All the results shown are divided by monopole contributions of the no-wiggle Kaiser model. 

total redshift power spectrum. These higher /i contributions have large amplitudes with differing signs 
[29], which would suggest that we might not be able to rely on our expansion in low-fc any longer, 
although FoG resummation can still help here. From figures we can also see relative contributions to 
the total monopole and quadrupole power from /i 2 and /x 4 terms. We see that at scales larger than 
k ~ 0.15h/Mpc /i 4 term contributes with 5-10% (depending on the redshift) to the total power of 
monopole and with 15-30% for the quadrupole. For the quadrupole all the remaining power comes 
form the /i 2 term while for the monopole the /x 2 term constitutes 25-35% of power and the rest 
comes from isotropic Poo term. To reduce the dynamical range we again divide monopole results by 
the no-wiggle monopole Kaiser term and the quadrupole results by the no-wiggle quadrupole Kaiser 
term. 

4 Conclusions 

In this paper we use the distribution function approach to redshift space distortions (RSD) that 
decomposes RSD into moments of distribution function. Our goal is to model the terms that contribute 
to the redshift space distortions using perturbation theory. We first repeat the derivations presented 
in [28] , explicitly deriving the decomposition of the moments into helicity eigenstates, based on their 
transformation properties under rotation around the direction of the Fourier mode. We give the 
explicit forms of correlators of the moments of distribution function and their angular dependencies. 
The leading order contributions to RSD can be classified in terms of their angular dependence, with 
the lowest order being /i 2 and /x 4 , where \i is the angle between the Fourier mode and the line of sight. 
There are three terms contributing to fj, 2 and seven terms contributing to /i 4 . We evaluate all of 
these terms using the lowest order PT (one loop) and compare them to simulation results. For some 
terms adding two loop contributions proves to be important and we extend our models to include the 
relevant contributions. Also, for some of these terms standard PT is not sufficient and we propose 
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Figure 13. Quadmpole moment Pi a is plotted at four redshifts, z = 0.0, 0.5, 1.0 and 2.0. To fj, 2 (red) term 
of P ss expansion we add /i 4 term (blue). We show PT (dashed lines) results, improved velocity dispersion 
model presented in this paper (solid lines), simulations up to fi 4 contributions (diamonds and stars), and 
reference simulation results (points). Resummed FoG model from equation 3.62 (black dot-dashed line) and 
linear Kaiser model (black dotted line) are also shown. All the results shown are divided by quadrupole 
contributions of the no-wiggle Kaiser model. 



physically motivated ansatz that goes beyond the loop analysis. These are based on the small scale 
induced velocity dispersion effects which multiply the long range correlations, such as density-density 
or density-velocity correlations. Such ansatz has a free parameter, small scale velocity dispersion, 
which cannot be modeled using PT. We found that a number of these terms have the same value, 
but also that not all velocity dispersions should be equal. We developed a halo model to describe the 
hierarchy of these terms and shown that the model can qualitatively explain the simulation results. 
Our analysis systematically accounts for all of the PT terms at one loop order and the small scale 
dispersion parameters, while necessary for a good description of RSD, have physically motivated 
values. In this sense our model goes beyond previous analyses [9, 10], which include some, but not 
all of the PT terms and which often treat FoG parameters as fitting parameters without a physical 
meaning. 

The dominant term to RSD is the /i 2 term and its dominant contribution is the momentum 
density correlated with the density. This term can be written in terms of a time derivative of the 
power spectrum [28] and so can be modeled using dark matter power spectrum emulators. Two other 
terms contribute to /j? , the vector part of the momentum density-momentum density correlation, and 
the scalar part of energy density-density correlation. We find that they affect RSD at a 10% level 
already at k ~ 0.05h/Mpc. The energy density-density correlation term is the dominant nonlinear 
effect, is negative for all scales and thus reduces the total /z 2 power. It is related to the Fingers-of-God 
effect. This term contains velocity dispersion term which cannot be modeled in PT and requires a 
free parameter in the model. 

The next angular term has /i 4 dependence and there are seven terms that contribute to the total 
power spectrum, of which one, scalar part of Pn, contains a linear contribution that does not vanish 
on large scales. We evaluate all of these terms in PT. Some of these terms are well modeled by PT, 
while others also require velocity dispersion type parameters. With these the modeling achieves some 
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level of precision compared to the simulations, but is still limited in the dynamic range, with an error 
of about 5% at k ~ 0.2h/Mpc at z = 1. 

Our ultimate goal is to develop accurate models of RSD that can be applied to observations. 
We observe galaxies, not dark matter, and understanding the physical processes that lead to RSD in 
dark matter is just the first step towards the goal of understanding the RSD in galaxies. The results 
presented here are only a rough guide for the challenges awaiting us when applying these techniques 
to the data, but there are some lessons learned that are likely to be valid also for galaxies. One is 
the importance of velocity dispersion effects, which dominate our model uncertainties on large scales. 
The good news may be that the velocity dispersion effects, which are the main source of the modeling 
difficulties in this paper, may be smaller for galaxies than for the dark matter, specially if a sample of 
central galaxies can be selected. We also expect that the halo model for computing velocity dispersion 
should be applicable to galaxies as well. However, galaxies also have additional challenges not present 
for dark matter: galaxy biasing will introduce additional scale dependent effects in redshift space that 
will need to be modeled, even if there are no such scale dependent biases in real space [28, 43]. The 
success of modeling the RSD and extracting the cosmological information from it depends on our 
ability to model these galaxy biasing and velocity dispersion terms. We plan to address some of these 
issues in the future work. 
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A Components of moments of distribution function 

Starting from equation 2.6 we can chose some basis, for example Cartesian, to express scalar product 
h ■ q. It follows 

T h L (x) - ^ f d 3 p/(x, P ) (h tPl /ma) L , (A.l) 
P J 

where summation over i = 1,2, 3 is implied. Using the multinomial theorem: 

{x 1 + X 2 + ... + X m ) n = klUt J l . .fc^ l 1 ^ 2 ---^' ( A ' 2 ) 

fci+fe+...+fc m =n 

it follows that 

If(x) = — J2 V^ur^ lh 2 h t l I dV(x,p)^ 1 p 2 fe2 ^/(ma) L . (A.3) 

P k 1+ k 2 +k 3 =L K l' K 2-'«3- J 

Neglecting velocity dispersion and anisotropic stress second rank tensor, and similar higher rank 
tensors contributions (cr l - J = 0, . . .) we are left with 

J <fpf{^)pfpfpf/{ma) L - (1 + *(x))«* 1 t;* 2 t;*», (A.4) 

where Vi is given in equation 2.4. Returning this back into equation A.3 and using multinomial 
theorem again we get 

l£(x) = (1 + <5(x)) k£^^ hlVl)kl{h2V2)k2 ^ V ^ = (1 + 5(x))(h ' V(X))L ' (A ' 5) 

ki + k 2 +ks=L 

Thus we have retrieved result of equation 2.7, and choosing h = f we get equation 2.13. 
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B Decomposition of in spherical tensors 

In this section we want to retrieve, starting from equation 2.15, equation 2.16. Let us consider the 
object T^(x) as defined in 2.6, which can actually be constructed by contracting all the components 
of rank L of tensor f __ A (equation 2.15) with unit h vectors. Fourier transforming this object 
gives us simply 

T-(k) = ^/^ /(k)P )f^ i (B.l) 



P 

Since we have translation symmetry it follows 

(x|x + r)=£(r) =>. (k |k') = (2ir) 3 P(k)5 D (k - k'), 

where it is implied that we take correlation of a general function of similar form like equation 2.6. This 
enables us to work with each Fourier mode separately, and add them appropriately in the end when 
we discuss the power spectrum. By symmetry of the problem we may choose a reference frame where 
z-axis is along h vector. Since spherical harmonics form a complete set of orthonormal functions and 
thus form an orthonormal basis of the Hilbert space of square-integrable functions, we can expanded 
/(k, p) in that frame as a linear combination, 

oo m—l 

f(k,p,8,<t>)=J2 E /r(k,p)3WM), (B.2) 

1=0 m=-l 

where p is the amplitude of momentum. Let us now consider the transformation properties of /, m (k, p) 
under rotation around the z-axis. We can think of rotation by some angle ip (i.e. <p' — cf> + ip) in two 
ways 

oo m=l 



/(k,p,M') = /(k,P,M + ^) = E E /T(M%».(M) 

1=0 m=-l 
m—l 



oo m=l 



1=0 ', 



From rotation properties of spherical harmonics it follows that /™(k,p) transform as 

/r(k, P )' = e "'^/ ; m (k,p), (b.3) 

so it is an eigenstate of an helicity operator id /dip, with helicity eigenvalue, or simply helicity, m. 
A quantity with helicity is called a scalar, that with helicity m = 1 is called a vector and that 
with m = 2 a tensor, but the expansion goes to arbitrary values of m. It is possible to do similar 
considerations for arbitrary rotation and so it can be shown that /™ transform as spherical tensors. 

In the chosen reference frame, using (h • p) = p z = p cos8, and inserting equation B.2 into 
equation B.l we obtain 

T h L (k) = E / dpp 2 fr(k,p)(p/am) L J dQ Y lm {9^)cos L 9 

l,m 

= Y^T l L > m (k)I lrn 5 m0 , (B.4) 

/ . rn 

where we have defined helicity eigenstates of moments of the distribution function 

T ; L < m (k) = 4tt^ J dp p 2 (p/am) L fr(Kp), (B.5) 
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and used abbreviation for the integral 



1 21 + 1 (I - m)\ 
47r [l + m)\ 



dx x L PJ n {x). (B.6) 



We have used definition of spherical harmonics Yi m {6, (j>) — J [;+")! Pi n ( c °sO)e m ^ , and x = cos(? 
abbreviation. Since in equation B.4 we have Kronecker delta 5 m o it is easy to evaluate integral 




2/ + 1 f 1 , r , . (-l) l L\ 121 + 1 f 1 , 

"ST ^ fl(af) = M^V "ST 

^^(l + C-l)*-') if I< i, 
if I > L, 



x L -\x 2 -l) 1 

(B.7) 



and here we have used nf = (7 ] ^-^7x77^ Jrr~^ ■ Collecting all that, equation B.4 becomes 

V t / (2 (-^+'+3)J 

^(k)= 2 J^nflf'°(k). (B.8) 

(l=L,L-2,...) 

Since we are working in z || h frame it is apparent that for some arbitrary k, the angular dependence 
is contained in T ; L '°(k) spherical tensors. The goal now is to disentangle the angular dependence 
from the radial. The procedure depends on whether one is using active or passive interpretation of 
rotation transformation. Let us first look at active interpretation. Then the completely contracted 
tensors, like the one we are dealing with equation B.8, can be obtained from the same one evaluated 
in k z = kz || h direction by rotating it in general k direction. Because T { L ' m are spherical tensors it 
follows 



(l=L,L-2,...) 



4tt 

-l 



= E E f^ntV«l{R)T^{kz), (B.9) 
(i=i,i-2,...) m=—l 

where T>(R) is the Wigner rotation matrix and "D^'mi^) ^ s ma trix elements. Using the well known 
relations, T>% m (Br l ) = V^ n ,(R) and Z><£>(0,0,O) = \/lfe*? ro (M), we get 



m—l 

lf(k) = J2 E nfT^ m {k)Y lm {6A\ (B.10) 

(l=L,L-2,...) m=—l 

where spherical harmonics now describe rotation form direction k back to h. On the other hand, 
using the passive interpretation we argue that (k) in frame z || h can be obtained by rotating it 
from z' || k frame, which is described by the same equations as before B.9. Note now that we were 
able to express result in terms of T ; z, ' m (fc), just a function of amplitude k, and all angular dependence 
is given with spherical harmonics which now describe the angular dependence in h direction seen 
from z' || k frame. Now setting simply h = f along a line of sight direction we get result 2.16, where 
cos# = r ■ k/fc = fi. 

Finally we show that in decomposition form B.10 we retrieve the same number of independent 
components (L + l)(L + 2)/2 as inferred from symmetries of equation 2.15, since for even L, i.e. I = 2n 
we have 

, , v^, N L #(# + 1) (L + l)(L + 2) 
^(21 + 1) = £(4n+ 1) = 1 + - + 4 2V2 2 = ^ >j >-, 

n=0 n=0 
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and for odd L, i.e. I = 2n + 1, we have 



(L-l)/2 



(£-l)/2 



L - 1 



L-l /-L-l 



£ (2/ + l) = J2 (4n + 3)=3 + 3— +4 



+ (i + l)(i + 2) 



n=0 



n=0 



C Conjugation properties of Pz,L'(k) 

In this section we investigate the conjugation properties of Pll 1 (k) functions. Starting from the 
condition that overdensity field <5(x) and velocity field v(x) are real valued fields, it follows that for 
Fourier space fields <5(k), 0(k) and fy(k) we have /*(k) = /(— k). This is valid also for more complex 
fields like 



Pn(k) 



d 3 qid 3 q 2 . ..d 3 q r , 



(2Tr) 3n 

If we compute conjugated field we get 
d 3 q 1 d 3 q 2 . . . d 3 q. 



/i(qi)/2(q 2 ) ■ • • /n(q„)^(k - q x - q x 



(C.l) 



Pn(k) 



(2ti-) 3 ™ 
d 3 qid 3 q 2 . . . d 3 q, 



(27T) 



3n 



A*(qi)/ 2 *(q 2 ) • • ■ /:(q„)^(k - qi - qi - . . . - q J 

/i(-qi)/2(-q 2 ) ■ • • /«(-q„)<5 £ '(k - q x - q x - . . . - q„) 



= p„(-k). 

From equation 2.7 it follows that T|* L (k) = lif (— k), so for correlator we have 

'T|f (k)7]f L '(k')) = (rf (-k)^ L (-k')) = (if (k)7|f L (k') 
thus we have Pi£/(k) = P£, L (k). So, for sum of two correlator we can write 

T,f(k)I]T z '(k / )) + (if' (k)n L (k')\ - 2Re(T,f (k)r,r L '(k') 



(C.2) 



(C.3) 



(C.4) 



D Integrals I(k) and J(k) 

Here we define integrals I nm (k) and J n m(k) used in previous chapters: 
I nm (k) = J J^jE /nm(k,q)Pz,(g)Pz,(|k- q|) and J„ m (fc) = 

where we define kernels / nm (k, q), and use r = q/k and a; = k • q/ (fcg): 



(2tt) 3 W UJ g 2 : 



(D.l) 



/oo(k, q) 

/io(k, q) 
/o 2 (k, q) 
/l 2 (k,q) 
/22(k, q) 
/3o(k, q) 
/i 3 (k,q) 
/ 32 (k,q) 



7:r+3r — lOrx 
14r(l+r 2 -2rx) / 

x(7x+3r-10ra; 2 ) 
14r 2 (l+r 2 -2r:E) ' 
(x 2 -l)(7x+3r-10rx 2 ) 

14r(l+r 2 -2rx) 2 ' 
(x 2 -l)(7x-r-6ra; 2 ) 

14r(l+r 2 -2rx) 2 ' 
x(7x— r— 6rx 2 ) 
14r-2(l+r 2 -2rx) ' 
1 — 3x 2 — 3rx+hrx 3, 

r 2 (l+r 2 -2rx) ' 
4rx+3x 2 — 6rx 3 — 1 

2r 2 (l+r 2 -2r:c) ' 
(l-z 2 )(2-12ra-3r 2 + 15r 2 3; 2 ) 
r 2 (l+r 2 -2rx) 2 ■ 



/oi(k,q) 

/n(k,q) 
/2o(k,q) 
/2i(k,q) 
/o 3 (k,q) 
/ 3 i(k,q) 
/ 23 (k,q) 
/ 33 (k,q) 



(7x+3r-Wrx 2 )(7x-r-e,rx 2 ) 
(14r(l+r 2 -2ra:)) 2 ' 

( 7x— r— 6rx 2 \ 
14r(l+r 2 ~2rx) J ' 
(2a+r-3rz 2 )(7a;+3r-10r3: 2 ) 

14r 2 (l+r 2 -2rx) 2 ' 
(2x-\-r—3rx 2 ^ (7x—r—6rx 2 ) 

14r 2 (1+r 2 -2rx) 2 ' 
(l-a: 2 )(3ra;-l) 
r 2 (l+r 2 -2r:r) ' 
(l-2rx)(l-x 2 ) 
2r 2 (i+ r 2 -2rx) ' 

3(l-ir 2 ) 2 
(l+r 2 -2ra;) 2 ' 

-4+12x 2 +24rx-40rx 3 +3r 2 -30r 2 x 2 +35r 2 x l 
r 2 (l+r 2 -2ra;) 2 
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Also we have kernels j nm (r): 
1 /12 



3oi(7 
SioO" 
ffn (r 
9Q2 (r 
9w(r 



3024 \r 2 
1 /24 



- 158 + 100r 2 - 42r 4 



(r 2 -l) 3 (7r 2 + 2) In 



3024 V r 2 
1 



202 + 56r 2 - 30r 4 + — (r 2 - l) 3 (5r 2 + 4) In 



' r + 1 

r + 1 " 
r- 1 



1008 
1 /12 



38 + 48r 2 - 18r 4 + -(r 2 - l) d In 
r 



r + 1 



|r-l| 



1008 V r 2 



82 + 4r 2 - 6r 4 + ^ (r 2 - l) 3 {r 2 + 2) In 



M^^ + l) (3^ -1^+3)-!^ -I) 4 In 



r + 1 
r + l 



r-1 



1 



9 



-2 (9 - 109r 2 + 63r 4 - 27r 6 ) + — (r 2 - l) (3r 2 + l) In 



672 \r 2 

and all the rest vanish in next to leading order regime. 



r + l 



(D.2) 
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